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ABSTRACT 


Equivalent impedance of a rough surface is determined by solving the problem of 
radiowave propagation over rough random surface by using the parabolic equation method. 
We consider horizontal polarization and assume the rough surface, defined by discrete 
points, to be perfectly conducting. This assumption has minimum effects for frequencies in 
the VHF band and above. The essential point of this thesis is to solve numerically a 
Volterra integral equation of the second kind for the surface current. Because the 
computation becomes ven intense for long ranges, we present a multiple sections solution, 
where the solution procedure is successively repeated for smaller ranges. The numerical 
results are compared with results available in the literature. Finally, we apply the technique 
to determine the equivalent impedance of a sinusoidal random rough surface by comparing 


the resulting normalized field to the field of a 2-ray model over a constant impedance plane. 
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I. INTRODUCTION 

Low-grazing angle radiowave propagation over random rough surface is an 
extremely important topic in many areas of interest to the Navy, such as sea-to-sea 
communications or detection of low-flying targets over the sea surface. Depending on 
the rough surface characteristics, the direct field and the field reflected and diffracted by 
the rough sities will interact in different ways, changing the field strength at a certain 
point of interest. | Numerous analytical methods are available for predicting 
electromagnetic wave propagation over such surfaces, such as physical optics, geometric 
optics, normal mode analysis and combination of the above. However, a complex surface 
profile complicates the application of these methods. The Parabolic Equation (PE) 
method, where a Helmholtz equation is approximated by a parabolic equation, has been 
extensively used for the determination of forward propagation because of its numerical 
efficiency. When making this approximation we are assuming that waves travel 
predominantly in the forward direction and that the back-scattered waves are neglected. 
This situation is met in many practical propagation problems. Naturally there are 
situations where it is not possible to apply such an approximation. For example, when 
there is a large obstacle behind a receiver, the back-scattered field cannot be neglected. 
The use of the parabolic partial differential equation method allows us to calculate the 
field at one location in terms of the field at a previous location, utilizing a marching 
scheme, resulting in a more efficient computational method. 

In this thesis we numerically solve a Volterra integral equation of the second kind 
to predict fields over a rough surface. We consider horizontal polarization and a 


perfectly conducting rough surface to simplify the formulation. The later assumption will 





have minimum effects for frequencies in the VHF band and above. Then we apply the 
technique developed to niin the equivalent impedance of sinusoidal rough surfaces. 
This is done by comparing the resulting normalized field to the field of a 2-ray model 
over a constant impedance plane. A random surface with a fixed spectrum is simulated 
by translating the sinusoidal surface of given peak-to-peak height and period in a random 
fashion. This is done by including a phase for the sinusoidal surface and varying it 


uniformly over [0,27]. Finally we show how to generate the rough sea surface according 


to the Pierson-Moskowitz spectrum. 

Chapter II describes the theoretical formulation, the numerical solution used to 
solve the integral equation and the methodology for the equivalent impedance calculation 
when using parabolic equation. Chapter III presents numerical results, comparing them 
with practical cases of radiowave propagation over terrain, and calculates the equivalent 
impedance for a sinusoidal surface. Chapter IV describes how to generate the sea surface 


according to the Pierson-Moskowitz spectrum. Recommendations and conclusions are 


presented in Chapter V. 











Il. THEORETICAL FORMULATION 

In this chapter we first present a Volterra type integral equation of the second kind 
for the determination of the surface current on the ‘oi surface based on the parabolic 
equation approximation. Because a straightforward solution of a Volterra integral 
equation for long ranges becomes computationally burdensome, we present an alternate 
technique whereby (i) the range is divided into several sections, (ii) the field is calculated 
at vertical lines at the junction of two sections, and (iii) the field on a vertical line is used | 
as a source for the iiiahalie of current on the next section. The rough surface is 
defined by discrete points and a linear interpolation between consecutive points is 
assumed. Figure 1 shows a possible irregular surface. The terrain is assumed to be a 
perfect electric conductor (conductivity — infinite). A 2-dimensional case is considered 


where the source, geometry and fields are assumed to be invariant in the y direction. 


Propagation is assumed to take place in the xz plane and is mathematically described by 
the standard parabolic equation given by 


. 2 
Ou _ ok (2.1) 
ox 2k, dz 
where ’u’ is the reduced field variable and is related to the electrical field ’U’ through 


U(x, z)=e**u(x,z), 


20 


k, a5 is the wave number and / is the free-space wavelength. An e 


time | 


dependence is assumed throughout. 





The following assumptions are made: (i) the waves are assumed to predominantly 
travel in the positive x-direction, and (ii) back scattering of waves is ignored. Given a 
source field, f (z), at x=0 one wants to calculate the field at x, where x is the 
horizontal distance from the source and z is the height with respect to some reference 
level. An integral representation, similar to Kirchoff’s integral form, can be derived 
based on integrating the parabolic equation in the domain formed by two vertical liries, 


the rough surface, and an infinite boundary as shown in Figure 1. 
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Figure 1: Perfectly Conducting Rough Surface 





The field at the point (x, z) is given by [Ref. 1]: 


L fx 
u(x,z)= | f(MGolx, 230, man-> J 





2 
(EG, [x, z:£,n = g(é)dé 


def 


= U(X, Z)—U, (x, Z) (2.2) 








where i=¥—1 and G,[x, z;€,7] is the Green’s function for the basic parabolic equation 


i i(z—1)° ko 
G,[x, 236, =a’ 22-6) 2.3 


For later use we will also be interested in the plane-wave spectral representation 


in free space: 


of the Green’s function obtained by decomposing the field in terms of plane waves 


travelling back and forth along the z -axis: 


G[x,2€.n\= HS) e $) | eile 


**0 dp (2.4) 
where H(e) is the Heaviside step function. 

The first integral of Equation 2.2 is performed in the vertical direction using the 
source field f(7). Its lower limit is g, and its upper limit is infinite. However, when a 
numerical computation is required, it will be truncated at some finite height. This 
integral corresponds to the field produced by the source f (7) in free space. 


The second integral of Equation 2.2 is performed in the horizontal direction along 


the rough surface 7=g(¢). This integral corresponds to the field reflected and 
diffracted by the rough surface. The term a6) corresponds to the vertical derivative 
n 


of the field on the surface and is related to the normal derivative on the rough surface by 


ou 1 ou 


ee (2.5) 
a J1+[2’(x) 


where g(x) is the surface profile, and the prime denotes differentiation with respect to 


the argument. Whenever we refer to the term current we will assume it to mean the 


3) 














vertical derivative of surface field. An integral equation fo 


differentiating Equation 2.2 with respect to z and tending z to the surface. Taking the 


proper limiting values one gets [Ref. 1]: 





2 4] rons [x, ¢(x):0, 71d + 
7] 


45 a me een ea ge sO) ox, 8(x)5$, a(S) de (2.6) 


where + denotes a principal value integral. 


Equation 2.6 is a Volterra integral equation of the second kind where the current 
at a given point will depend only on the current at the previous points. Thus we notice 
the causal property of the Volterra integral equation. In order to solve it we need to use a 
marching procedure, where we calculate the current at a point based on the current at all 
other previous points. However, this current calculation becomes computationally very 
intense for large distances. In order to alleviate this problem, we present the multiple 
sections solution described on the next section. Once we have the current on the surface, 
given by Equation 2.6, we can find the field at any point (x, z) using Equation 2.2. 

A. MULTIPLE SECTIONS SOLUTION 

We can reduce the computation time by dividing the surface into a number of 
sections and using Equation 2.6 over the sub-intervals. Thus, instead of calculating the 
current on the entire surface between a transmitter and a location of interest, we can 


perform current calculations successively on shorter sections by making use of field data 


available on vertical lines. 











Given the field on a vertical line, the following procedure is performed in order to 
calculate the field on the next vertical line: 

Step 1: Calculate the incident current induced on the rough surface by the initial 
field. The incident current on the surface is given by the first term on the right hand side 
in Equation 2.6. 

Step 2: Calculate the total current on the rough surface. The total current on the 
surface is given by both terms on the right hand side of Equation 2.6. 

Step 3: Calculate the field at the next vertical line, utilizing Equation 2.2. 

We notice a trade off when applying this split-section method, illustrated on 
Figure 2. On one hand we reduce the calculation time needed to solve the Volterra 
integral equation since we are dealing with smaller distances in each individual section. 
On the other hand, it will be necessary to calculate the field on vertical lines at 
intermediate sections. There will be an optimum number of atihi where the total - 


computation time will be minimized. 
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Figure 2: Representation of the Multiple Sections Solution 





B. INITIAL FIELD SPECIFICATION 


In order to be able to band limit the initial data as well as to account for finite 


antenna beamwidths we use a Gaussian spatial distribution given by: . 


_ (=z)? 


= 20," 
u(O,z) = Ae (2.7) 
where z, is the transmitting antenna height, o, is it standard deviation for the Gaussian 


distribution and is related to the antenna beamwidth in the vertical plane through 


BW,,, =2tan™ al 


09, 
and A is the amplitude. Obviously, choosing a small o, (for example 4/2) makes the 
source field very concentrated around z,. Notice Figure 37 (upper figure) where we have 
the source field as a function of o,. 
Using this source field as f (7) in the first integral in Equation 2.2 we find that the 


free space field, which corresponds to the first integral of Equation 2.2, is given by: 


__(2=2,¥° 


U per qpaae(%sZ) =m (i (2.8) 


Given the source field described by Equation 2.8 we can also calculate its 


derivative at (x, z) with respect to z: 


(z—z,) 
ou (x, z) _ Ao, (z me: ) é ote) | (2.9) 
Oz 














C. NUMERICAL SOLUTION PROCEDURE 
1. Calculation of the Incident Current on the Rough Surface 
a. Calculation for the First Section 
The incident current on the surface is given by the first term on the right | 


hand side in Equation 2.6. 
a, A) : . 
Tine (x)= 2 | Ff in) Gol 8(2):0.mldn (2.10) 
80 | 


For the first section we have an analytical expression for f(7), given by 
Equation 2.7. Substituting z= g(x) , taking the derivative with respect to 77 in Equation 


2.3, and using this result and Equation 2.7 into the right hand side of Equation 2.10, we 
find the expression for the incident current on the rough surface. This expression is valid 


only for the first section where an analytical expression for the initial field is available: 


as 2 2| o2+— 
7,,,(x)=—2AGlee)- a] 2d, QL) 


b. Calculation for the Remainder Sections 
For the remainder steps we need to perform a numerical evaluation of the 
incident current, given by Equation 2.10. The use of trapezoidal rule is not suitable 


because the Green’s function term inside the integral in Equation 2.10 has a phase term of 
] a ae ’ 

the type —, and for small distances x it introduces numerical errors. As a result, a more 
x 


sophisticated integration technique is required. 





Making use of the spectral representation in Equation 2.4 we rewrite the 


expression for J, (x) as: 


imc 


ip*x 


J ige(X) = mt) [ F(n)[ petlrMe dp (2.12) 





Li) roy 
Interchanging the orders of integration and substituting 77 = 85 +& we get: 


Jipe(2)=—H(x) | peP*Or#'e dp] f(go+S)e ds 
0 


—oo 


def i 


< pts | 
—H(x) [ pe®lsrele % 3(p)dp (2.13) 


where 3(p)= | f(g, +é)eP' dE is the Fourier transform of the initial field. 


0 


Since the parabolic equation has a limited angular response [Ref. 2], we 


need to truncate the integration limits in Equation 2.13 at + p,,,,, Where P,,, 18 Of the 


order of kysin(6 nox ) where @_,, =15°. Because we are truncating the integral in 
Equation 2.13, we need to introduce a filter, such as a Hanning filter, in order to smoothly 
roll-off the function. The expression for J,,, then becomes: 
ip?x 


a ides _ip*x 
Tie(x)=— [W(p)pS(pye * ele #ldp (2.14) 


~ Pmax 


q 


where w( p) is the Hanning filter. The above integral will be evaluated numerically by 
the trapezoidal rule. The Hanning filter is an even function and its discrete form is given 


by: 


10 











47m 


W(p=mAp)= ; 
: on S| —msp<sM 
N ]8 


where M is the number of points used in the discrete Fourier transform. The final 


numerical expression for J,,. is 


Pax 
= 2 Yw(r,) (Pn )PnS(Pq Je 1 eivale hee aaa 


m=] 


where p,, =mAp. 


2. Calculation of the Total Current on the Rough Surface 
Once we calculate the incident current on the rough surface, the next step is to 


calculate the total current on the surface. The total current on the surface is represented 





0 
by Equation 2.6. We can say for the clarity of explanation that J(x)= e 
| | Z 





J(€)= (x;€) = 2-20 vee we can rewrite Eananon 2.6 as: 
x 
I(x) = Jig (X) + Fe f(g) SO mH ae (2.16) 
in. ri , HE 


where J, (x) is given by Equation 2.10. 


The second term on the right-hand side in Equation 2.16 can be broken in N sub- 


integrals: 


1] 





I(x) =J,,,(x +d a fa J(é) 222 8alx%g) c) greed) 3p 


n=l TE (“1)5 x-¢ 
=o (x)+¥ “oy (2.17) 
INC = Ti n 


no Ko ee 2 x: 
where I, = | J (€) "= £7050) o) im 5° E)ga ( UE. 
(n-]l)6 yx —€ 


Notice from Figure 3 below that the following notation has been adopted: 


n = 1,2,3,...,N +1 
x, =(n-1)d 
Xyay =NO=x 


x-d=x, =(N-lod 


g, = 8(x,) 

J, =J(x,) 
n 1 2 n n+l N+1 
x O 5 (n-1)5 nd Nd 


Figure 3: Nomenclature for the discretization process 


With reference to the Figure 4, we can approximate g,(x;¢)= so) as 
X= 


g(x) — 8(%,,.4) 
£4, =—————~.. where x ayy, =(n-26 and s(t,.y)=~ 


eV, 











approximate the variable £ in g,(x;&), located between x, =(n—1)6 and x,,, =nd by 


Xa =n, and the function value g(&), located between g, and g,,, by g(x y) 
nt Yo n+ fo 


2N+1 
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Figure 4: Approximation for g,(x;é) 
Letting € =76+(n—1)6, in Equation 2.17, where 0<7<1 in the n™ segment, 
and approximating J(é) linearly between J, and J,,, as J(Q)=(J,u4-J, \r+J,, we 
simplify the expression for /,, as: 


kd 
~i-—8 ant 


i052 (N+1-n) | e 
1, =¢,,v6e6? J.-J, e+ 5, —————— dt 2.18 
n B an € iI n+] \r (N+1—n)-T ( ) 


The result may be expressed in terms of the Fresnel integral as [Ref. 1]: 


‘in aE ihe {v- 1 \v+1—n)+J, a) Net )-Ft, I 


ik, O8 in 


Jri-J,) fe[ ee | 
_Ya=In), 2’ | the 28 —t, |p for 209, and J, =0 for g,, =0 (2.19) 
ik, Og a, 


13. 











+1x>0 








sen(x) =40,x=0 
—-1,x<0 
k,O(N—n 
ty = - San 

M4 

k, O(N +1—n) 

t, a | dn 
1 
and 
F(t)= fer dx is the complex Fresnel integral. (2.20) 


0 


Once we have the expression for J,,we can calculate the total current on the 


surface using Equation 2.17: 


2 Sn iT 4 a J, ) | 
J yar =F incon F | S88 an > Jiu wd, \N +1- n)+ J, an [F(t, )- F(t, + 
n=l - 0S dn 


where 


E: 
Ty =2sgn(8yu — 8 ),{7—X 
Ky 
m4} Jy { #3 Flt | 
«Feats Jo tae 2 i} (2.22) 
int; int, 2 | 


In Equation 2.22 we notice that 7,, =O whenever gy = 8yu1- 








Combining Equations 2.21 and 2.22 we have the final expression for the current 


on the surface: 
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a 1 | =3 Flt 
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3. Calculation of the Field 


Once we have the total current on the rough surface we can calculate the field at 


any x. In Equation 2.2 we notice that u,(x, z) corresponds to the incident field and 


u., (x, Z) corresponds to the field reflected and diffracted by the rough surface 


a. Calculation of the Incident Field 


The incident field, corresponding to the first integral on the right hand side 
of Equation 2.2 and given by u,(x,z)= [ f (MG, [x, z;0,7]d7 will be calculated using 
go 


the same technique as used for the calculation of J inc (X ). Itis given by 
px 


1 7 eee 
us(x,z)=— | S(p)e 2ko gir(z-8o) gy 





(2.24) 
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As before, a Hanning filter is introduced into the expression when the 
limits of the integral are truncated. The expression for u, (x, z) then becomes: 


Pmax px *x 


(x, z)=— —— | wl )S(p)e * ee 80) dp (2.25) 


ae 
Finally we can apply the trapezoidal rule to Equation 2.25, obtaining the final numerical 
expression for u, (x, z): 
Dax 


u, (x, z) z)= -~ yw p, B (p, Je 2H eiale-so) (2.26) 


n=l 
b. Calculation of the Field Reflected and Diffracted by the Rough 


Surface 


The field reflected and diffracted by the rough surface, corresponding to 


the second integral on the right hand side of Equation 2.2 and given by 


Ma(2)= hos [x E.n=a(é)laé (2.27) 


cannot be calculated by applying the trapezoidal rule directly because the Green's 


function has a singularity at x=¢. This singularity must be extracted apriori before 


performing the integration. We can rewrite the above expression for u, (x,z) as 





j x J(E) spose t 

2(x-€ 
u,{x,z)= —=-—¢ d 2.28 
, (x, z) ik Je ; (2.28) 


We replace [z—g(é)| by [z—9(x)+ 2(x)- g(E)| in Equation 2.28 and 


define: 


ik [2-2 oy! kole(x)-8 (EP 
K(x;€)=J(E)e (8) 5 rs) (2.29) 














so that 
K(x, x)= J(x)e™* Oe 6) (2.30) 


If g(x) is continuous K(x;&) is well behaved as x >. We can rewrite 











Equation 2.28 as: 
u,(x,z)= aa, x 
[z-e(x)P 
ik 2 
x” 2x-é) )—K(x: x) ik asad} 
x| K(x; x){ dé + K(n¢)— K5x) a2) dé (2.31) 
eter - a 


By achange of variables, the first term within square brackets of Equation 
2.31 can be solved in terms of the Fresnel integral resulting in [Ref. 1]: 


jkolz-s(z)P z—g(x)]’ 


ecskne Kae) Ks), ae a (2.32) 





u,(x,z)= ee 
0 


where: 


kolz-g(x)P 
1, (x, z) = 2Vxe 2% = 


—(1-i)/7k, aC 0-9F| a) 0.33) 


AX 


The integral in Equation 2.32 can now be evaluated by applying the 
trapezoidal rule to the second term on the right isd side, since the singularity at x =€ 
has been extracted. 

4. Fresnel Integral Calculation 
For small arguments (smaller than or equal than 0.1) the Fresnel integral can be | 


approximated by 
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F(x)=x+— ; (2.34) 


This expression is obtained by taking the limit when t->0 in Equation 2.20. Using this , 
approximation will reduce CPU time somewhat. For arguments greater than 0.1, the 


Fresnel integral will be calculated by a rational approximation given in [Ref. 3]: 


F(x)=C(x)+iS(x)=—(1+i)-[irl)+ se (2.35) 
where 
14+0.926x 
x)= —_ OO EX 
(x) 2+1.792x +3.104x? (x) 
[| a 
2+4.142x+3.492x? +6.670x° 
where le(x < 2x10~ 


D. DETERMINATION OF THE EQUIVALENT IMPEDAN CE 
Once we have determined the field on a vertical line distant x from the source, 


we can find the equivalent ground impedance that would produce the same field when 


solving the propagation problem using a two-ray model. 





For horizontal polarization, the reflection coefficient for plane waves incident 





from air on to a non-magnetic lossy, planar interface having an impedance 77 = 710 iS: 
€ TC 
sin —.j/e | 
[= e = (2.36) 


~ siny+Je,, 


where y is the grazing angle such that 


H.+H 


ar al ul 
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=sinWfor H,+H, <£ 


and 7, is the intrinsic impedance of the free-space. 
When considering the two-ray model the resultant field will be given by: 


kof; 
e 01 








Kor 
pe uqlebres (2.37) 
I. 


h ) 2 
where r, is the line-of-sight path, r, is the reflected path and u, is the free-space field. 


For large distances (d >> H, +H,) we can approximate r, ~ 7, for the amplitude 





term, and use 7, —7, = a for the phase term. We can rewrite Equation 2.37 as 


=14Teltol2-1) (2.38) 


0 


and express the reflection coefficient as: 


, i2koH,H, | | 
T= 7 _ i} 4 (2.39) 
Uy 


If we now combine Equations 2.36 and 2.39 we have the final expression for 


€,- as a function of the geometry and the normalized field: 


i2koH,H, 


Uu 

e ¢ +1-— 

le =; Ho 
€,. = Sm ch (2.40) 

e 4 —l+— 


Uy 


Figure 5 below illustrates the geometry mentioned above. 
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‘Figure 5: Reflection of Plane wave over Flat Surface 
Note that the heights H, and H, are relative to the equivalent flat plane, which iS 
normally taken to coincide with the mean of the rough surface. If only one section is 
used in the PE modeling we can take advantage of the Gaussian representation of the 
initial field. A somewhat better approximation to Equation 2.37 will then be: 


ey: 
Hott foe 
= x 






(2.41) 


j2koH Hr OT 
x ( ; (se =) 














Ii. RESULTS 

Having introduced the theoretical candienon in Chapter II, we will present here 
numerical results obtained with the algorithm developed and compare those with other 
techniques for some practical cases of radiowave propagation over terrain. We will be 
comparing the propagation factor (PF) (excess loss ue the free space case) obtained 
with our approach versus the propagation factor obtained from results available in the 
literature. A positive (negative) value of propagation factor implies gain (loss) with 
respect to propagation in free space. 
A. PROPAGATION OVER A PERFECT ELECTRIC CONDUCTOR PLANE 

Consider a transmitter 5 m above flat perfect electric conductor (PEC) ground 


(z,=5). One wants to calculate the PF on a vertical line distant 400 m from the 


transmitter, as shown in Figure 6. The maximum receiver height to be evaluated is 100 


1 


m. 





Figure 6: Propagation over a Flat PEC Surface 


The PF is defined as: 
PF [dBk 10log,, ula, z) (3.1) 
Uy (x, £ ) 
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where u(x,z) and u,(x, z) are defined in Equations 2.2 and 2.8 respectively. 
For a point source above a PEC plane an analytical solution can be obtained by 


image theory. Using the expression valid in the far field for x >> z,,z, we get 











Figure 7 shows a comparison between the two methods for a frequency of 
operation of 1 GHz and performing the computation using one section. For this range of 
400 m good agreement is seen for heights up to 60 m. This result is expected, since the 
2-ray model and the PE model use approximations that are valid for horizontal distance 
large compared to antenna heights. If we move the receiver vertical line further to the 
right on Figure 6, thereby, increasing the distance from the source, we will have better 
agreement at higher heights. 

On Table 3.1 we have the relationship between the CPU time (Pentium MMX 200 
MHz) and the number of sections used in the calculation of the field over a total distance 
of 400 m. The code was implemented in MATLAB. The column % — refers to a 
relative measure to the maximum measured time. We notice that as we change the 
number of sections from one to two and from two to three the CPU time decreases. 
However, as we — from three to four sections the CPU time increases. Of course 
changing the section length will make the CPU time change too. Compare in Table 3.1 
cases two and three, both using 2 sections; sages four and five, both using 3 sections and 


cases six and seven, both using 4 sections. Figure 8 illustrates this dependence. These 


numbers are influenced by many factors, some dependent on the set-up of the problem 

















(number of points on the vertical direction, for example) and some dependent on the 
computer configuration (available memory, for example). 

We 4is6 notice that it is not possible to choose the sections too small, since it 
would nei that waves would have to travel at large angles to reach large vertical heights 


and the parabolic equation has an angular limitation around 15°. Referring to Figure 9 


we notice that the results for the cases 3, 5 and 7 are smoother than 2, 4 and 6 


respectively. This is because the fields on a vertical line get more accurate as the 
horizontal distance from the source increases. 

We can make a better check to our numerical solution by comparing it with the 
solution for Gaussian sources instead of the point sources used above. For a Gaussian 


source, described by Equation 2.7, an analytical expression for the total field over a 


ground plane is: 


where the first term on the right hand side corresponds to the free-space field and the 
second term corresponds to the field generated by the image source. The PF will be given 


by 


(3.2) 





Figure 10 shows current and the propagation factor for Case 2 (d1 = d2 = 200 m). 


There we notice a small oscillatory behavior for the current. Compare this figure with 
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Figure 7, where the calculations were performed in one single section. It is seen that the 
results are smoother. Below we offer an analysis of this behavior. 

Figure 11 compares the propagation factor at the end of the first section in Case 2 
(d1 = d2 = 200 m), when using our numerical solution and using the image theory with 
Gaussian source as described by Equation 3.2. We notice such a good agreement, even 
for large heights, that it is hard to distinguish the two curves. Figure 12 compares the 
magnitude and the phase of the field at the end of the first section. For the sake of clarity 
we present Figure 13 that contains basically the same information as Figure 12 but over a 
smaller vertical scale. There we compare our numerical solution with the solution 
obtained using the image theory described above. We notice a very good agreement, 
which indicates that the field at the end of the first section is being calculated correctly. 
As a result, we conclude that although the multiple section solution improves the speed of 
the calculation, it introduces some numerical errors, since the PF for Case 1 and Case 2 
are not the same, and the only difference between the two cases is that in Case 2 we 


performed the computation in two sections. In Section 3.C we will see that as we 


decrease the length of the first section the results become worse. 
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Figure 7: Current (Upper Figure) and Propagation Factor (Lower Figure). Transmitter 


antenna height is 5m. Flat earth, frequency = 1 GHz, A, = - 
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Figure 8: CPU Time as a Function of Number of Sections for the PEC Flat Surface Case 
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Figure 9: Propagation Factor as a Function of the Number of Sections (x-axis: PF in dB, 


y-axis: Receiver Height in meters). 
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Figure 10: Current (Upper Figure) and Propagation Factor (Lower Figure) for the Case 2 


of Flat Surface. Transmitter antenna height is 5m. Frequency = 1 GHz, A, = a 
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Figure 11: Propagation Factor, Calculated at the End of the First Section of Case 2 of Flat 


Surface. 
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Figure 12: Magnitude and Phase of the Field at the End of the First Section for Case 2 of 


Flat Surface. 
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Figure 13: Magnitude and Phase of the Field at the End of the First Section for Case 2 of 


Flat Surface. 
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B. PROPAGATION OVER PEC KNIFE-EDGE 


Consider an ideal absorbing knife-edge of height h, =5 meters on a conducting 


plane, as shown in Figure 14, where we have a = b = c = d = 5O m. The transmitter 


antenna height is 3 meters (z, =3) and it is located at point A. One wants to calculate 


the PF at a vertical line located at B, distant 200 m from the transmitter antenna. 


receiver 
height 





Zt 


Figure 14: Perfectly Absorber Knife-Edge Between the Transmitter at A and a Vertical 
Line at B, Where the Propagation Factor is Being Evaluated. Both are over Perfectly 


Conducting Ground. 


We compare our ne solution with a four-ray model of knife-edge 
diffraction theory [Ref. 4]. In the parabolic equation computations we replace the knife- 
edge with a triangular hill having a base of 100 m (b = c = 50 m in Figure 14) and a peak 
height of 5m. The numerical results are plotted in Figure 15. We notice that the total 
current and the incident current are the same for the first 50 m (corresponding to portion 
a), since we have flat surface for this part of the terrain. We notice a good agreement for 


low heights, since the four-ray model assumes large distances. If we now move both the 


triangular peak and the vertical line B to the right, making a = d = 150 m and b=c = 50 
32 














m = Figure 14 we have the numerical results shown in Figure 16. Moving the triangular 
peak and the vertical line B further to the right, making a = d = 250 m and b=c =50m 
on Figure 14, we have the results shown on Figure 17. Comparing Figure 15, Figure 16 
and Figure 17 we notice progressively better agreement with increasing heights, since 


both models compared provide more accurate results for longer horizontal distances. 
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Figure 15: Total Current and Incident Current (Upper Figure) and Propagation Factor 
(Lower figure) for the Triangular Hill with a = b = c = d= 50 m, Peak Height =5 m and 


Antenna Height = 3 m. Frequency = 1 GHz. Calculation in one section. A, = A. 
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Figure 16: Total Current and Incident Current (Upper Figure) and Propagation Factor 
(Lower figure) for the Triangular Hill with a= d = 150 m, b = c = 50 m, Peak Height = 5 


m and Antenna Height = 3 m. Frequency = 1 GHz. Calculation in one section. A, =A. 
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Figure 17: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower figure) for the Triangular Hill with a = d = 250 m, b = c = 50 m, Peak Height = 5 


m and Antenna Height = 3 m. Frequency = 1 GHz. Calculation in one section. A, = A. 











C. PROPAGATION OVER A PEC SINUSOIDAL SURFACE AND EQUIVALENT 
IMPEDANCE DETERMINATION 

Once we tested our solution for two important test cases, we can use the model for 
a rough surface. The only restriction we enforce is on the slope of the surface: as it 
increases, the waves are forced to travel in higher angles and the parabolic equation has 
an angular limitation of approximated 15° [Ref. 2]. In this section we will analyze the 
case of a sinusoidal surface. Initially we will evaluate the propagation factor and verify 
how fast we can march, 1.e., how the results are sensitive to changes of the step size. 
Finally, we will determine the equivalent impedance for low grazing angles. 

1. Propagation over a PEC sinusoidal surface 


Consider a sinusoidal surface defined by 


s()=“e{1-coo +0 (3.3) 


where A,, is the peak to peak amplitude, x is the horizontal distance from the source, 
T is the period, @ an arbitrary phase and g(x) the amplitude on that particular horizontal 


distance. 

In Table 3.2 we have the relationship between the CPU time (Penttum MMX 266 
MHz, MATLAB code) and the number of sections used in the calculation of the field over 
a total distance of 800 m on the sinusoidal surface with the following parameters: 


A, =1m, T= 30 m, @=0. The distances di, d2, d3 and so on are defined in Figure 18. _ 
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Table 3.2: CPU time for Propagation over PEC Sinusoidal Surface 





Figure 19 ae the incident current, the total current and the propagation factor 
when performing the calculation in one section. Figure 20, Figure 21 and Figure 22 show 
the incident current, the total current and the propagation factor when performing the 
calculation in two sections, gradually decreasing the length of the first section (d1) and 
increasing the length of the second section (d2). We notice that as we decrease the length 
of the first section and increase the length of the second section we have a disagreement 
among the results. Comparing cases 2 and 3 with case 1 the disagreement is not 


significant. However, case 4 and case 1 have completely different results. 
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Figure 23 shows the incident current, the total current and the propagation factor 
in three sections. Figure 24 shows the incident current, the total current and the 
propagation factor in nine sections. 

Comparing Figure 19 and Figure 24, where we have the calculation in one section 
and nine sections respectively, we notice that the currents and the PF become very noisy 
for the last case. However, we have the computation time reduced to 41 % of the original 
one. In contrast to the case of flat surface, we observe that for propagation over rough 
surface we need to have many small sections to reach the point such that the 
computational time becomes minimum. This does not happen for the flat surface case 


because we had all g,,’s equal to zero and we did not have to calculate the Fresnel 


integral when evaluating the J, ’s (refer to Equation 2.19). 
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Figure 19: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for the Case 1 of Sinusoidal Surface. Antenna Height = 5 m. Frequency 


= 1 GHz. 
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Figure 20: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for the Case 2 of Sinusoidal Surface. Antenna Height = 5 m. Frequency 


= 1 GHz. 
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Figure 21: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for the Case 3 of Sinusoidal Surface. Antenna Height = 5 m. Frequency 


= 1 GHz. 
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Figure 22: Total Current and Incident Current (Upper Figure) and Propagation Factor 





(Lower Figure) for the Case 4 of Sinusoidal Surface. Antenna Height = 5 m. Frequency 


= 1 GHz. 
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Figure 23: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for the Case 5 of Sinusoidal Surface. Antenna Height = 5 m. Frequency 


= 1 GHz 
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Figure 24: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for the Case 10 of Sinusoidal Surface. Antenna Height = 5 m. Frequency 


= 1 GHz 
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2. Sensitivity to the step size, A,, and to the standard deviation,o, , of the 


Gaussian source field 


Consider propagation over a sinusoidal surface with the following parameters: 
z,=5m, A, =1m, T=30m, x=300m, 0=0° and frequency of 1 GHz (A =0.3m). 


Figure 25, Figure 26, Figure 27, Figure 28, Figure 29 and Figure 30 shows the incident 
current, the total current and the PF for A, equal to “ (0.15m), A (0.3m), 24 (0.6m), 


4A (1.2m), 8A (2.4m) and 16A (4.8 m) respectively and the source field with o, =0.1. 


Figure 31, Figure 32, Figure 33, Figure 34, Figure 35 and Figure 36 shows the incident 


current, the total current and the PF for A, equal to “ (0.15m), 2 (0.3m), 24 (0.6m), 


4A (1.2m), 84 (2.4m) and 164 (4.8m) respectively and the source field with 


We notice that for o, =0.1 we have good results for both the total current and 
the PF for A, up to 2A. For o, =0.4 we have good results for the total current for 
A up to 16/, although even for small A, (A/2) we have problems with the PF for high 
heights. This happens because the free-space field becomes very small for those heights 
due i the Gaussian distribution and numerical errors are introduced. In Figure 37 we 
notice that for o, =0.1m the free space field has almost constant magnitude at the 
receiver location, while for o, =0.4m the field magnitude decays as we go higher. The 
field becomes too engi for high heights and numerical errors are introduced when 


computing the PF. 
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Figure 25: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for Propagation over Sinusoidal Surface. A, = oe 0,= Et 
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Figure 26: Total Current and Incident Current (Upper Figure) and Propagation Factor 
A 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =A, 0, = ; 
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Figure 27: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =21, 0, = ; 
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Figure 28: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =4/1, 0, = - 
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Figure 29: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =8/1, 0, 3 
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Figure 30: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =161, o, = : 
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Figure 31: Total Current and Incident Current (Upper Figure) and Propagation Factor 


(Lower Figure) for Propagation over Sinusoidal Surface. A.= “ 70. = “ 
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Figure 32: Total Current and Incident Current (Upper Figure) and Propagation Factor 


Piss 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =/, 0, 3 
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Figure 33: Total Current and Incident Current (Upper Figure) and Propagation Factor 


=s4 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =2A, Oo, : 
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Figure 34: Total Current and Incident Current (Upper Figure) and Propagation Factor 


ele 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =4/, 0, : 
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Figure 35: T otal Current and Incident Current (Upper Figure) and Propagation Factor 


see 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =8/, 0, : 
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Figure 36: Total Current and Incident Current (Upper Figure) and Propagation Factor 


ae 


(Lower Figure) for Propagation over Sinusoidal Surface. A, =161, 0, 3 
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Figure 37: Normalized Source Field (Upper Figure) and Free Space Field (Lower Figure) 


at x = 300 m. Transmitter Antenna Height = 5 m. Frequency = 1 GHz 











3. Equivalent Impedance Determination 

Once we have calculated the propagation factor for the sinusoidal surface, we can 
find the equivalent impedance of the sinusoidal surface for low grazing angles. If we 
compare Figure 7 and Figure 19 we notice that for small receiver heights (= low grazing 
angles) both curves have approximately the same shape. This can give us an intuitive 
explanation that under these circumstances we can represent the PEC sinusoidal surface 
by the lossy flat surface. 


Figure 38 shows the real part and the imaginary part of the normalized equivalent 





impedance for the sinusoidal surface using the following parameters: z, =1.5m, 


S 


re 


A, =Ilm, 7=30m and x=300m. We ran the problem 50 times with @ uniformly 


distributed in the interval [0,27]. We can see that as the grazing angle approaches zero, 
the real part tends to 1.69 and the imaginary part tends to 4.58. Notice the dashed line in 
the first plot on Figure 38, which represents the best linear fit when minimizing the sum 
of the squared error. If we now mn the flat surface problem with these values we obtain 
the same results for PF as the sinusoidal surface when averaged over many realizations. 
Notice the bottom plot on Figure 38. Figure 39 siowe results when running 100 
realizations. We notice that as the grazing angle approaches zero, the real part tends to 
2.16 and the imaginary part tends to 4.56. 

In Figure 40 we ran 50 realizations with the parameters described above 


(z,=1.5m, A, =1m, T=30m, x=300m and @ uniformly distributed in the interval 


{0,27]), but we changed the maximum receiver height to 78 m, or, equivalently, the 


maximum grazing angle to 15°. Then, using the equivalent impedance found when 
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performing the fitting for low grazing angles (real part tending to 1.69 and imaginary part 
tending to 4.58) we calculated the equivalent PF for propagation over flat lossy surface. 
In the upper figure we notice a linear behavior for both the real aT imaginary part of the 
equivalent impedance for a maximum grazing angle of approximately 7°. In the lower 
figure we notice a good agreement between the propagation factor when performing 50 
realizations and the equivalent flat lossy surface for a maximum height of approximately 
15 m. This result is expected, since as we observed in Figure 7 and Figure 19, we have 
approximately the same shape only for small receiver heights. In other words, the 
substitution of the sinusoidal PEC surface for a lossy flat surface, when comparing the 


PF, is valid only for low grazing angles. 


For extremely small grazing angles, the errors in = are amplified by Equation 
Uo | 


. This is the cause of the large swings in Figure 38, 





2.39, due to the presence of 


Figure 39, and Figure 40 (upper figures) at low grazing angles and is also the cause of 
discrepancy seen in the PF in the PF in Figure 38, Figure 39, and Figure 40 (lower 


figures). Notice that we have smaller swings for one hundred realizations than fifty 


realizations. 
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Figure 38: Upper Figure: Equivalent Impedance of Sinusoidal Surface with Period = 30 

m and Peak to Peak Amplitude = 1 m. Total Distance = 300 m. Frequency = 1 GHz. 

Antenna Height = 1.5m. Ax =0.3m. Fifty realizations. Lower Figure: PF for the PEC - 


Sinusoidal Surface Above and for the Equivalent Flat Lossy Surface. 
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Figure 39: Upper Figure: Equivalent Impedance of Sinusoidal Surface with Period = 30 
m and Peak to Peak Amplitude = 1 m. Total Distance = 300 m. Frequency = 1 GHz. 


Antenna Height = 1.5m. Ax = 0.3m. One hundred realizations. Lower Figure: PF for 


the PEC Sinusoidal Surface Above and for the Equivalent Flat Lossy Surface. 
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Figure 40: Upper Figure: Equivalent Impedance of Sinusoidal Surface with Period = 30 
m and Peak to Peak Amplitude = 1 m. Total Distance = 300 m. Frequency = 1 GHz. 
Antenna Height = 1.5m. Ax = 0.3m. Fifty realizations. Lower Figure: PF for the PEC 


Sinusoidal Surface Above and for the Equivalent Flat Lossy Surface. 
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IV. SEA SURFACE GENERATION 
The sea spectrum can be modeled in several different ways depending on the 
parameters that one tries to represent. In this chapter we will use the Pierson-Moskowitz 


spectrum, which represents the sea spectrum as a function of the wind speed and is given 


by [Ref. 5]: 


(4.1) 


where @=8.1x10°, 8 =0.74, g =9.81 m/s” , U =wind speed (m/s) at a height of 19.5 
m, k= = is the wave number and A is the surface wavelength. 


Figure 41 shows the spectrum for wind speed of 5m/s and 6m/s. 
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Figure 41: Pierson-Moskowitz Spectrum for Wind Speed of 5m/s and 6m/s 
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Realizations of the random surface will be generated at points 


x, =nAx (n=1,...,N) according to [Ref. 6]: 


vA 
it : IK ;x 
fle, )=— SF(K,)e™”, | (4.2) 
L oN 
where, for j20, 
[N(0,1)+ iN (0,1)] N 
: 
F(K,)=,[2aLW{K,) (4.3) 
_.N 
N(0,1), j= 0, 


and, for j<0, F(K,)=F'(K,). In Equation 4.3 x=, where L is the surface 


length and each time N (0,1) appears it represents an independent sample taken from a 
zero mean, unit variance Gaussian distribution. Equation 4.2 is computed with a fast 
Fourier transform. 

Figure 42 (upper figure) shows a possible surface realization for a wind speed of 
5m/s. Since the parabolic equation has a limited angular response [Ref. 2] we can avoid 
high slopes smoothing the surface, i.e., changing a local point by the average between 
itself and some neighbors points. Figure 42 (lower figure) shows the realization of the 
upper figure when replacing any point by the moving average of itself and four 
neighbors, two to the right and two to the left. When making use of this smoothing 
process we practically do not change the statistics of the surface. Figure 43 (upper 


figure) shows a comparison between the inverse Fourier transform of the noisy spectrum 


and the autocorrelation of the rough surface for the wind speed of 5m/s. Figure 43 














(lower figure) shows the same comparison, but replacing the autocorrelation of the rough 
surface by the autocorrelation of the smoothed rough surface. 

Defining the correlation length ’]’ as the length such that 0.5x/ equals the area 
under the autocorrelation function (only for positive lags) and d,as the dian 
between two consecutive points on the rough surface, we need to have 

L>>I and I>>d, 
in order to have a truly uncorrelated surface. In the example above, where we had the | 
_ wind speed of Sm/s, we had L= 300m, /=21.2 m and d, =0.29 m. 

Figure 44 shows a surface realization (top figure) and the respectively smoothed 
surface realization (bottom figure) for a wind speed of 20m/s. Notice that the surface 
heights increase substantially when compared to Figure 42. Figure 45 (upper figure) 
shows a comparison between the inverse Fourier transform of the noisy spectrum and the 
autocorrelation of the rough surface for the wind speed of 20m/s. Figure 45 (lower 
figure) shows the same comparison, but replacing the autocorrelation of the rough surface 
by the autocorrelation of the smoothed rough surface. For this case we of wind speed of 


20m/s we have L= 2000 m, 1 = 469.9 m and d, =3.9m. 
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Figure 42: Sea Surface Realization for Wind Speed of 5 m/s (Upper Figure) and 


Smoothed Sea Surface (Averaging Five Points) for Wind Speed of 5m/s (Lower Figure). . 
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Figure 43: Normalized Fourier Transform of the Noisy Spectrum and Rough Surface 
Autocorrelation (Upper Figure) and Normalized Fourier Transform of the Noisy 


Spectrum and Smoothed Rough Surface Autocorrelation (Lower Figure). Wind Speed = 


5m/s 
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Figure 44: Sea Surface Realization for Wind Speed of 20m/s (Upper Figure) and 


Smoothed Sea Surface (Averaging Five Points) for Wind Speed of 20m/s (Lower 


Figure). Length = 4000m 
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Figure 45: Normalized Fourier Transform of the Noisy Spectrum and Rough Surface 
Autocorrelation (Upper Figure) and Normalized Fourier Transform of the Noisy 
Spectrum and Smoothed Rough Surface Autocorrelation (Lower Figure). Wind Speed = 


20m/s 
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Vv. CONCLUSIONS AND RECOMMENDATIONS 

In this thesis, a numerical method to model propagation over random rough 
surface using the parabolic equation was implemented and tested. The parabolic 
equation, where the Helmholtz equation is approximated by a parabolic equation, is a 
quasi-full wave method and takes into account aspects such as forward reflection, 
diffraction, refraction and surface wave propagation. However, it ignores back- 
seatioang: Because of this assumption, it is possible to use a marching scheme, where 
the field at a certain point is calculated uniquely based on the field at a previous point, 
resulting in a very efficient computational method. Besides, this method seems to be the 
only practical method for predicting radiowave propagation over long ranges (thousands | 
of wavelengths) and over a wideband of frequencies. 

Since our technique uses a Volterra integral equation of the second kind, the 
numerical computations become very intense for large distances. In order to minimize 
the computation time, we introduced a multiple sections solution, where the numerical 
solution scheme is successively repeated for saiailes ranges. Given the field on a vertical 
line between two sections we calculate the current on the surface between these two 
vertical lines and then calculate the field on the next vertical line. When using this 
multiple sections solution we substantially decrease the computation time, although we 
noticed that when sections become too smal! numerical errors were introduced. 

We tested our numerical techniques for some test cases encountered in the 
literature such as propagation over a PEC flat surface and over an ideal absorbing knife- 


edge. The agreement between the numerical solution and these test cases was good. 
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Next we tested our numerical solution for propagation over a sinusoidal surface. 
We noticed that for low grazing angles the propagation factor has a shape similar to the 
one when propagation takes place over a flat lossy surface. ‘This fact served as 
motivation for modeling the difficult problem of propagation over a rough surface at low 
grazing angles by propagation over an equivalent impedance plane. We simulated a 
random surface with a fixed pera by translating the sinusoidal surface of given peak 


to peak height and period in a random fashion, with a uniform distribution in [0,27]. 


Then we compared the resulting normalized field to the field of a 2-ray model over a 
constant impedance plane. For low grazing angles the agreement was excellent. Finally 
we showed how to generate rough sea surfaces according to the Pierson-Moskowitz 
spectrum, which represents the sea spectrum as a function of the wind speed. 

This thesis represents an effort to predict radiowave propagation over a random 
rough surface such as the sea. The next step should be to predict the propagation factor 
for propagation over the sea surface. Finally, because the solution of the Volterra integral 
equation takes significant amount of CPU time during recursions, the source code should 
be implemented in a programming language that deals with for loops efficiently, such as 
C or FORTRAN. In this thesis we used MATLAB, which is now to be very efficient for 


matrix manipulation, but very inefficient for loops. 
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